function [StdMean, StdCovar] = ecmnstd(Data, Mean, Covar, Method)
%ECMNSTD Estimate standard errors for mean and covariance of incomplete data
%	NUMSERIES x 1 column vector of standard errors of estimates for each
%	element of a mean vector Mean and NUMSERIES x NUMSERIES matrix of standard
%	errors of estimates for each element of a covariance matrix Covar, where
%	Mean and Covar are maximum likelihood estimates derived from Data, which
%	has NUMSAMPLES independent identically-distributed samples of NUMSERIES
%	random variables assumed to be multivariate normal with missing data
%	indicated with NaNs.
%
%	StdMean = ecmnstd(Data, Mean, Covar);
%
%	[StdMean, StdCovar] = ecmnstd(Data, Mean, Covar);
%	[StdMean, StdCovar] = ecmnstd(Data, Mean, Covar, Method);
%
% Input:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a
%		NUMSERIES-dimensional random vector. Missing values are indicated
%		by NaNs.
%	Mean - NUMSERIES x 1 column vector of maximum likelihood parameter estimates
%		for the mean of Data.
%	Covar - NUMSERIES x NUMSERIES matrix of maximum likelihood covariance
%		estimates for the covariance of Data.
%
% Optional Inputs:
%	Method - String to indicate method of estimation for standard error
%		calculations. The default method is 'hessian'. The methods are:
%		'hessian' - (default) Use the Hessian of the observed negative
%			log-likelihood function. This method is recommended since the
%			resultant standard errors incorporate the increase uncertainties
%			due to missing data.
%		'fisher' - Use the Fisher information matrix.
%
% Output -
%	StdMean - NumSeries x 1 column vector of standard errors of estimates for
%		each element of the mean vector Mean.
%	StdCovar - NumSeries x NumSeries matrix of standard errors of estimates for
%		each element of the covariance matrix Covar.
%
% WARNING: If calculating standard errors associated with the covariance matrix
%	Covar, i.e., two output arguments, then this routine is VERY slow.
%
%	See also: ecmnfish, ecmnhess, ecmnmle, ecmnobj 

%	Author(s): R.Taylor, 4-21-2005
%	Copyright 2005 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $   $Date: 2005/06/17 20:23:25 $

% Step 1 - check arguments

if nargin < 3
	error('Finance:ecmnstd:MissingInputArg', ...
		'Required inputs Data, Mean, and Covar are missing.');
else
	[NumSamples, NumSeries] = size(Data);

	Mean = Mean(:);
	if ~all(size(Mean) == [NumSeries, 1])
		error('Finance:ecmnstd:IncompatibleMean', ...
			'The mean vector Mean has wrong dimensions.');
	end
	if ~all(size(Covar) == [NumSeries, NumSeries])
		error('Finance:ecmnstd:IncompatibleCovar', ...
			'The covariance matrix Covar has wrong dimensions.');
	end
end
if nargin < 4
	Method = 'HESSIAN';
else
	Method = upper(Method);
	if ~any(strcmp(Method,{'HESSIAN','FISHER'}))
		warning('Finance:ecmnstd:UnknownMethodString', ...
			'Unknown Method string. Will use default HESSIAN.');
		Method = 'HESSIAN';
	end
end

if nargout < 2
	MatrixFormat = 'MEANONLY';
else
	MatrixFormat = 'FULL';
end

% Step 2 - initialization

if any(sum(isnan(Data),1) == NumSamples)
	error('Finance:ecmnstd:TooManyNaNs', ...
		'One or more data series has all NaN values.');
end

if sum(sum(isinf(Data)))
	error('Finance:ecmnstd:InfiniteValue', ...
		'One or more infinite values found in data.');
end

[CholCovar, Flag] = chol(Covar);

if Flag > 0
	error('Finance:ecmnstd:NonPosDefCovar', ...
		'The covariance matrix is not positive-definite.');
else
	CholCovar = inv(CholCovar);
	InvCovar = CholCovar * CholCovar';
end

% Step 3 - calculate Fisher information or Hessian and invert

if strcmp(Method,'FISHER')
	if strcmp(MatrixFormat,'MEANONLY')
		ParamCovar = Covar;
		Count = sum(~all(isnan(Data),2));
		ParamCovar = (1.0/Count) .* ParamCovar;
		StdMean = sqrt(diag(ParamCovar));
		return
	else
		ParamCovar = ecmnfish(Data,Covar,InvCovar,MatrixFormat);
	end
else
	ParamCovar = ecmnhess(Data,Covar,InvCovar,MatrixFormat);
end

Count = sum(~all(isnan(Data),2));

MeanCovar = inv(ParamCovar(1:NumSeries,1:NumSeries));
MeanCovar = (1.0/Count) .* MeanCovar;

if strcmp(MatrixFormat,'FULL')
	CovarCovar = inv(ParamCovar(1 + NumSeries:end,1 + NumSeries:end));
	CovarCovar = (1.0/Count) .* CovarCovar;
end

% Step 4 - pull standard errors out of diagonal elements

StdMean = sqrt(diag(MeanCovar));

if strcmp(MatrixFormat,'FULL')
	StdCovar = zeros(NumSeries,NumSeries);
	k = 0;
	for i = 1:NumSeries
		for j = 1:i
			k = k + 1;
			StdCovar(i,j) = sqrt(CovarCovar(k,k));
			StdCovar(j,i) = StdCovar(i,j);
		end
	end
end
